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An experimentally validated approach is described for fast axisymmetric Stirling engine 
simulations. These simulations include the entire displacer interior and demonstrate it is 
possible to model a complete engine cycle in less than an hour. The focus of this effort 
was to demonstrate it is possible to produce useful Stirling engine performance results in 
a time-frame short enough to impact design decisions. The combination of utilizing the 
latest 64-bit Opteron computer processors, fiber-optical Myrinet communications, dynamic 
meshing, and across zone partitioning has enabled solution times at least 240 times faster 
than previous attempts at simulating the axisymmetric Stirling engine. A comparison of 
the multidimensional results, calibrated one-dimensional results, and known experimental 
results is shown. 

This preliminary comparison demonstrates that axisymmetric simulations can be very 
accurate, but more work remains to improve the simulations through such means as mod- 
ifying the thermal equilibrium regenerator models, adding fluid-structure interactions, in- 
cluding radiation effects, and incorporating mechanodynamics. 

Nomenclature 


7 Porosity of the Medium 

r e ff Stoke’ s Tensor 

Pf Fluid Density 

p s Solid Density 

v Fluid Velocity Vector 

Ef Fluid Energy 

E s Solid Energy 

h Sensible Enthalpy 

J Diffusion Flux of Helium 

k e ff Effective Thermoconductivity 

p Pressure 

Sj Fluid Enthalpy Source Term 

T Temperature 


I. Introduction 

Great strides have been made in the last decade in expanding the designer’s analysis options 
for building improved free-piston Stirling engines 1 ranging from thermal and structural dynamic 
analysis, 2 ’ 3 Computer-assisted-design (CAD) software for complete three-dimensional part 
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design, 4 one-dimensional thermodynamic cycle analysis 5,6,7,8,9,10,11,12,13,14,15,16,17,18, and optimization, 
and one-dimensional system mechanodynamic simulation. ’’’’’’ A review of the progress 
in Stirling analysis shows a steady relaxation of assumptions/limitations imposed on the analysis 
over the past century. 26 

Initial two-dimensional computational fluid dynamic (CFD) simulation work began over 
twenty years ago and focused on specific regions of the engine but were not time-efficient and the 
projects were put on hold. 6,27,28,29,30 And other researchers attempted to build specific in-house 

Q 1 rsrs rsA q r 

research codes with mixed success at modeling specific components. ’ ’ ’ ’ ’ 

More recently, multidimensional multiphysics simulations that allow for realistic Stirling 
engine system solutions now exist and include the following capabilities that are potentially 
useful for Stirling convertor design: 37,38,39,40 


• Structural 

— linear or nonlinear static or transient deformation 

— modal, spectrum, harmonic, and random vibration dynamics 

— linear or nonlinear buckling 

• Thermal 

— steady- state, transient, conduction, convection, and radiation 

• CFD 


— Steady-state/transient, compressible/incompressible, laminar /turbulent, forced/natural convec- 
tion, conjugate and radiation heat transfer, laminar to turbulent transition modeling, scalable 
parallel performance 

• Electromagnetics 

— electrostatics, magnetostatics, low-frequency electromagnetics, current conduction, circuit analy- 
sis, harmonic, transient, modal, scattering, perfect electric and magnetic conductors, impedence 
boundaries, near and far electromagnetic field extension, frequency selective surface, antenna 
radiation patterns, and specific absorption rate. 

• Coupled Physics 

— acoustics-structural, electric- magnetic, fluid-structural, fluid-thermal, elect romagnetic-thermal, 
piezoelectric, piezoresitive, thermal-electric, thermal-structural, electromagnetic-thermal-structural, 
and eletrostatic-structural 

Recently, attempts at utilizing these commercially available options have focused on only the thermal 
dynamic cycle simulation because it is difficult to achieve a useful solution in a reasonable time, let alone 
attempt a more comprehensive simulation. For example, previous attempts at whole engine simulation 
with these multiphysics commercial codes ended in unconverged solutions . 31,41 Although a recent simplified 
axisymmetric solution 75 shows reasonable energy balance convergence, the results are not validated yet by 
experiment. 

The first reported successful three-dimensional whole Stirling engine simulation was reported by Mahkamov 42, 43,44 
in the United Kingdom, but due to export control regulations, limited information and verfication is avail- 
able. Apparently a gamma-type Stirling engine was simulated using commercially available software with 
good experimental agreement. It is not clear how much detail was simulated nor how practical the results 
are to engine design. The second reported successful three-dimensional Stirling engine simulation was by 
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Zhang 45 in China in which it took approximately 3 months to reach a steady-harmonic converged solution 
on a simplified, academic geometry that did not include the appendix gap, solid walls, or internal displacer. 

A multi-dimensional analysis has an important place in Stirling engine design as is discussed at length 
in Ref. 46 and outlined below: 

• verify the one-dimensional results, 

• properly simulate inherently three-dimensional turbulence-including transition, 

• provide empirical heat transfer and friction factor coefficients for complex geometries before the hard- 
ware is built 

• integrate all the parts into a CAD system and test for structural and relative motion issues in the 
overall design 

• assist experiments with determining hard to reach data 

• provide a fluid-structure interaction capability 

• generate linear reduced order models for controller development 

• solve high-power Stirling applications in which the one-dimensional flow assumption breaks down 

• identify areas of excessive flow losses due to unintended dead zones, recirculation zones, dissipative 
turbulence and other losses such as: 47, 48,49 

1. Inefficient heat exchange and pressure loss in the regenerator, heater and cooler, 

2. Gas spring and working space loss due to hysterisis and turbulence, 

3. Appendix gap losses due to pumping and shuttle effects, 

4. Mixing gas losses from nonuniform temperature and flow distributions perpendicular to primary 
engine flow axis, 

5. Conduction losses from the hot to cold regions 

6. Losses due to combined radiation, conduction and convection in hollow spaces 

7. And in general, inaccurate loss representations due to use of 1-D flow design codes to account 
for flow and heat transfer through area changes (between components) where phenomena such as 
flow separation and jetting from tube into regenerator may occur. 

The current multiphysics simulation tools are not optimized for Stirling analysis as discussed in Ref. 50 
For example, recently developed oscillating time advance strategies 51,52,53,54,55,56,57 and high order spatial 
differencing could make modeling oscillating turbulent transition a reality. But given the cost of developing 
new codes, it is desirable to fully utilize current technology first. 

This paper will present the first successful U.S. axisymmetric Stirling engine simulation incorporating an 
actual engine’s geometry including appendix gap, displacer seals (inner and outer), solid walls, and displacer 
interior. Only the flexures, piston seals, bounce space, radiation shields, and heat exchanger fins were not 
included for this initial work, but reasonable agreement with experiment has been achieved despite these 
omissions. 


II. Description of the Problem 

The dual opposed configuration shown in figure 1 58 A9,60,6i j g p e i n g developed for multimission use 
(i.e., for use in atmospheres and space), including providing electric power for potential missions such as 
unmanned Mars rovers and deep space missions. 62,63 
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Only the Stirling engine part of the convertor 
(Fig. 2) is simulated multi-dimensionally although 
one could anticipate the entire convertor may one 
day be simulated using modified multiphysics soft- 
ware. 

An example of the desired full Stirling engine 
simulation is shown in Fig. 3. It is anticipated that 
full 3D simulations will provide a level of geometric 
and flow detail necessary for further design improve- 
ments. For example, hardware experiments have shown that large performance gains can be made by varying 
manifolds and heat exchanger designs to improve flow distributions in the heat exchangers. 59,64,65,66 


Figure 1. Dual Opposed Stirling Convertors Reduce 
Vibration (Schreiber) 



Figure 2. Stirling Engine Part of Convertor Simulated Only 


The kinds of flow and geometry that occur in a 
Stirling engine for which a modeling technique must 
address are as follows: 49,48 

1. Oscillating flow which changes the effective 
conduction, flow friction coefficients and heat 
diffusion. 67 

2. Low mach number flow (no shocks), 

3. Compressible flow due to enclosed varying vol- 
umes and heat transfer effects, 

4. Laminar, Transitional, and Turbulent flow 
with Reynolds numbers from 100 to 10000 
(based on various length scales pertinent to 
the flow region), 

5. Conjugate heat transfer, thermal dispersion 
and local thermal nonequilibrium in the 
porous media regenerator, 



Figure 3. Example Transient Simulation of Stirling 
Convertor-Colored by Temperature 


6. Micron to Millimeter scale geometry 6 
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7. Sliding Interfaces in the appendix gap 

8. Deforming Flow Regions Due to Compression 
and Expansion 

III. Current Multi-Dimensional Stirling Cycle Analysis Tools 

At the level of multidimensional analysis, relatively little has been completed because the third order 
analysis, including Sage, 5 and to some extent, DeltaE, 13 and MS*2 Stirling Cycle Code 11 are faster and for 
the most part have been adequate engineering design tools. However, to improve efficiency, reliability, and 
supplement experimental testing (to understand and reduce losses, identify component temperatures, assist 
in sensor calibration, correlation, and placement) it will likely require a better understanding of the actual 
flow features and heat transfer throughout the engine. Some of the primary multidimensional analysis tools 
are shown below. 

1. Modified Computer Aided Simulation of Turbulence (CAST) 

The modified CAST code 69 is based upon the Semi-Implicit Method for Pressure-Linked Equations (SIM- 
PLE) 70 method but is restricted to two dimensions. 71,72,73 It was modified to include oscillatory boundary 
conditions and conjugate heat transfer but does not allow for complicated geometry nor sliding interfaces. 
It has been used to model Stirling components but has not been extended to a whole engine simulation tool. 

A pressure-splitting technique was added 74 to reduce the computational requirements. It was based 
on separating the thermodynamic and hydrodynamic pressures so that these widely varying scales could 
be solved with less round-off error and better efficiency. And the non-acoustic form of the equations was 
later incorporated into the code since the fast acoustic waves and small dimensions of the engine essentially 
resulted in ” incompressible” like flow behaviour and removing the acoustical pressure disturbance could 
enhance computational speed and accuracy. 32 

2. CFD-ACE 

This commercial code has been used to model a two-dimensional representative Stirling engine 31,75,33,45 and 
it has been used to model an actual engine. 41 Full validation of these simulations has not been completed. 

It is also based upon the SIMPLE technique. The regenerator is not currently modeled correctly since 
thermal equilibrium is assumed between the gas and solid, but as shown later, this error may not be large 
compared to the overall engine performance. The code does support sliding interfaces on parallel computers 
in three-dimensions and the company is introducing sliding, overlapping Chimera grids which should enable 
better appendix gap modeling. Moreover, the highly compressed volumes occuring at both ends of the 
displacer may be more easily modeled with Chimera grids. This finite volume code can utilize both structured 
and unstructured grids. 

3. Fluent 

This commercial code is based upon the SIMPLE and PISO methods. It currently has similar regenerator 
modeling limitations as CFD-ACE. It does however have a sliding interface for two-dimensional axisymmetric 
flows that can be used for appendix gap modeling on parallel computers and will be discussed in more detail 
later. It is also finite volume based and can utilize both structured and unstructured grids. 

It has been used with mixed success by industry for Stirling engine simulations 76 by several commercial 
manufacturers. It was also used by Mahkamov 43 in what appears to be the first full three-dimensional 
simulation, but details are embargoed. Also, in the related business of cryocooler modeling it has been used 
with reasonable success on simplified problems 77, 78,79,81,82,83,84,85,86,87,88 And this report has utilized it to 
provide a detailed axisymmetric simulation. 
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I STAR-CD 

The Simulation of Turbulent Flow in Arbitrary Regions (STAR) 38 code also uses SIMPLE and PISO methods. 
Its companion product (STAR-HPC) is the parallel computer version. It also has sliding interfaces and 
deforming mesh capability. It has been used in the related field of internal combustion engine piston modeling, 
and some Stirling engines and regenerators have been modeled with it. 89,90 

5. CFX 

This code also uses SIMPLE and PISO methods on unstructured grids. It also has sliding interfaces imple- 
mented, but no Stirling engine modeling with this software has been publicly published. 37 It currently seems 
to have the greatest potential for full multiphysics Stirling convertor simulations because of it’s fluid-structure 
interaction, magnetic flux analysis, and structural/thermal stress analysis capabilities. It’s dynamic meshing 
capabilities are not ideal for highly compressed volumes as occur in the Stirling engine since only the spring 
analogy mesh adaption is currently available. 

6. Others 

While there are other research codes, they are usually limited to modeling only specific regions of the Stirling 
engine such as the regenerator, or the displacer. 80, 81, 15,91, 92 

IV. Regenerator Modeling 

A very important specific area of modeling difficulty is the regenerator. Since the regenerator (depending 
upon one’s definition of effective) has roughly 3 to 40 times more effective heat transfer than the heater, 93 
any inefficiency of the regenerator represents a significant loss for the entire Stirling engine. Hence, any 
numerical losses or inaccuracies in this region will disproportionately influence the entire Stirling simulation. 

Also, the geometrical shape of the matrix elements in the regenerator is important and as shown in 
Figs. 4(b-f), many shapes are possible. For practical calculations a Representative Elementary Volume 
(REV) is often used as shown in Fig. 4(a). An REV requires some closure assumptions that can be partially 
derived from experiment and potentially from microscale (pore level) numerical simulations. 

Since the size and shape of the regenerator fibers do affect Stirling engine performance, numerically 
simulating a regenerator has been pursued at both microscopic (direct) and macroscopic (REV) levels of 
fidelity. 6,94,95,96,97,98,99, 100 An example simulation utilizing Fluent by Harvey from Georgia Tech is shown 
in Fig. 5. 

However, our initial approach at whole engine simulation simply adds source terms to the governing 
equations and does not take into account the different temperatures of the solid matrix and the fluid. The 
source terms serve to enforce a pressure drop according to the Darcy-Forcheimer 67 equation. And only a 
single energy equation is used, although it does include an effective solid/fluid average conductivity 68 and 
an effective solid/fluid average energy. 

The continuity equation is not modified, but the following source term is added to the momentum 
equations (for i = x,y, orz ): 


Si = - (j^ v i + C 2 ^PVmagVi'j 


( 1 ) 


where permeability K = 4.08e 10 m 2 , C 2 = ^ '= and inertial coefficient, Cf = 0.178. These coefficients 
were determined experimentally 101 and used without modification here. And the fluid/solid averaged energy 
equation can be written as: 
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(a) Representative Elementary (b) Random 

Volume 


(c) Wire Mesh 



(d) Perforated Disk 


(e) Foam Metal 


(f) Sintered Glass (Failed) 






Figure 4. Regenerator Geometry 
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( 2 ) 


V. Dangers of Component Modeling 

Modeling the entire engine, as demonstrated herein, is preferable to modeling one component at a time due 
to the oscillatory nature of the flow. A component requires known unsteady boundary conditions for an entire 
cycle. Since a given boundary may have flow going in both directions, the direction of the characteristics 
at these boundaries are difficult to define. Incorrect specifications of boundary conditions result in solution 
instability, nonconvergence of solutions, and/or convergence to inaccurate results. It is well known that even 
for steady harmonic flows, non-physical reflections can occur at inflow and outflow locations. 102 In essence, 
by artificially truncating the domain to model a single engine component an artificial boundary condition 
is required. However, the exact specification of artificial/non-reflecting boundary conditions in nonlinear 
turbulent flow is an important and as yet unresolved area of research. 103 

For example, a subsonic inflow has four characteristics entering and hence four physical boundary con- 
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(a) Schematic 


(b) Pressure 


(c) Pathlines 


Figure 5. Flow Through Wire Mesh Screen Regenerator 


ditions need to be specified, typically total pressure, total temperature, and two flow angles, allowing only 
one of the velocity components to be free (which is extrapolated from the interior flow solution). A sub- 
sonic outflow has only one characteristic entering and hence only one physical boundary condition should be 
specified, for example, static pressure. 

Simultaneous inflow /out flow along each engine component’s interface is clearly demonstrated by exam- 
ining flow velocities near zero (from -1 m/s to 1 ms). Regions of blue and red in Fig. 6 show borders of 
mixed inflow/outflow in the axial direction. For example, the piston compression space in Fig. 6(a), the 
seal exit in Fig. 6(b), the cooler exit in Fig. 6(c), and the appendix gap entrance in Fig. refhg:oscillating(d) 
all show mixed inflow/outflow borders at this moment in time. Moreover, these mixed boundaries moving 
during the cycle, further complicating the numerical boundary treatment. If a component boundary crosses 
any of these indicated areas, the boundary conditions will not be properly defined. Of course in the case of 
commercial software users the code will attempt to stabilize the solution with essentially guessed boundary 
conditions, but the solutions will be inherently non-physical. 

Repeating the same velocity test, but this time examining velocity close to zero in the y-direction, we 
see a new set of boundary lines that will pose difficulties resulting in non-physical component solutions in 
Fig. 7. 

A. Adjust Components within a Full Engine Simulation 

Instead of pursuing the mathematically ill-posed problem of component modeling in complicated oscillating 
flow, it is easier, more efficient, and more reliable to determine a whole engine steady- harmonic solution 
with the procedures described herein. The engine components can be easily examined and it only takes a 
few hours to test the effect of component changes since the overall engine has the correct steady-harmonic 
temperature distribution. 

At the same time, it is possible to numerically determine the one-dimensional equivalent friction factors 
and heat transfer terms in the new components required by one-dimensional analysis tools (e.g. Sage). In 
this way, a one-dimensional analysis be made more accurate while not requiring the ill-posed boundary con- 
ditions that multidimensional component modeling entails. 46 Finally, it is worth noting that one-dimensional 
analysis does not have difficulty with component modeling because a boundary is always either an inflow or 
an outflow, never both. 
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(a) Entire Engine 


(b) Seal Region 




(c) Compression Space 


(d) Heater Manifold and Appendix Gap 


Figure 6. 


Oscillating Axial-Direction 
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(a) Entire Engine 


(b) Heater Region 



(c) Appendix Gap Region 


(d) Seal Region 


Figure 7. 


Oscillating Radial-Direction 
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VI. Approach for Fast Simulations 


Achieving practical computer simulations of entire Stirling engines depends upon judicial use of computer 
hardware, software engineering, grid generation, and algorithm design. A combination of all four areas is 
presented next that enables axisymmetric whole engine simulations to be accurately completed in less than 
a week. 


A. Computer Design 

An important area new to Stirling analysis is the use of parallel computers for faster simulations. NASA 
Glenn was one of the first organizations to pioneer computer clusters back in 1992. Many universities have 
assembled their own low cost computers, 104 and many states have supercomputer consortiums 105 to assist 
in designing and operating small computer clusters that utilize commercial off-the-shelf (COTS) hardware 
and software for low cost systems. These systems are referred to as Linux Beowulf clusters 106, 107 and now 
represent the majority of the top 500 supercomputer systems in the world. 108 Recently, many commercial or- 
ganizations have begun assembling ” Beowulf’ -like clusters including Dell, 109 Angstrom, 110 LinuxNetworx, 111 
Appro, 112 and Microway 113 as shown in Fig. 8. They cost a bit more, but provide comprehensive health 
monitoring and controls that significantly reduce the adminstration and integration costs. The latest systems 
can house 260 processors in a single 7 foot tall tower. 



(a) Angstrom Cluster 


Rear View 



1 master node 


SdaeeTar 
2 switches 


Power reed 


(b) Appro Cluster 



(c) Linux Networx Cluster 


Figure 8. Appro, Angstrom, and Linux Networx Clusters 

Our sixteen dual-Opteron (32) processor Microway cluster is shown in Fig. 9. This was selected because 
last year it represented the best value for performance and it integrated well with the commercial software, 
Fluent. The Opteron chip, at the time, was the only chip that could run 32-bit and 64-bit applications with 
an exceptional price/performance ratio. All the whole Stirling engine results shown here were calculated 
using it. 

A very important consideration when designing your computer cluster is the communication bandwidth 
and latency since commercial fluid simulation software utilizes older numerical approaches that require a large 
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(a) Front View 


(b) Dual CPU Board 


Figure 9. Microway Cluster 


number of frequent communications between the compute nodes. In Fig. 10, two of the best communications 
network interface card (NIC) and network switches are shown. The Infiniband system can perform 10 Gb/s 
bi-directional communication rates or about 20 times faster than the Myrinet in our cluster. Our simulation 
utilized Myrinet communications because last year Infiniband was not a proven technology. 

B. Software Decision 

We chose Fluent as our computational platform because, at the time, it was the only commercial vendor 
supplying an axisymmetric sliding interface for parallel computers. Presently, as mentioned in Sec. Ill, 
the other vendors are now also providing similar capability. It is important to always have at least two 
commercial solvers available because of the invariably large number of bugs in their software. Generally, a 
bug is more easily found if the same simulation is repeated with two different codes. For example, Fluent 
version 6.1.22 properly handles heated moving solids, but the new release (version 6.2.16) did not as identified 
by Ibrahim. 114 The latest version 6.2.19 has been corrected. 

Another reason for utilizing at least two commercial codes is the lack of control the user has in the 
solvers utilized. Generally, the commercial codes all utilize slight variants of older, well established lower 
order numerical techniques. The slight variants are the only independent way to examine the effect of a 
given solver and to look for possible defects. Of course developing one’s own software offers more freedom, 
but this can be risky in the short-term. 

Our second platform utilizes CFD-ACE and reasonable solutions have been achieved with it. It’s primary 
limitation has been the need to work in three-dimensions when solving cases including axisymmetry sliding 
interfaces. Another strength of this software is it includes a wide range of multiphysics capabilities. 

A third very interesting platform is CFX because of it’s integrated multiphysics capability enhanced by 
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(a) Myrinet Card 


(b) Myrinet Switch 


(c) Infiniband Card 


(d) Infiniband 144 port 
Switch 


Figure 10. Communications 


the merger between ANSYS and CFX. The only apparent limitation of the software is grid layering is not 
yet available. Grid layering is defined by the automatic addition or removal of a row of cells as the mesh 
domain shrinks or expands and will be discussed later. 

Other packages such as STAR-HPC are very capable as well. However, the commercial cost of these 
packages are very expensive and it is not practical to test all of them. 

Finally, full three-dimensional simulations are still quite challenging and the parallal scalability of the 
commercial software is somewhat limited due to the implicit nature of their steady flow solvers. Better 
approaches are currently available 50 but have not yet been applied to the Stirling cycle analysis. 

C. Grid Design 

When using commercial codes, it is primarily the grid quality that determines the speed and accuracy 
of the solution. Ideally, one could simply read in a three-dimensional CAD geometry and a grid would 
automatically be created. That is in fact the implied goal of computational fluid dynamics code developers. 
But for complicated geometry as occurs inside a Stirling engine with all it’s moving parts, the technology 
currently available still requires a great deal of human intervention. Taking several weeks to develop a grid 
is not unusual and is often required to get a reasonable solution. For example, the rate at which conjugate 
heat transfer can be solved is directly related to the smoothness and quality of the grid at the interface. This 
is particularly true in Stirling engines when the heat capacity of the solid is vastly diffent than the Helium. 
If short-cuts are taken in developing the grid, the price is often vastly longer solution times, if a solution is 
produced at all. 

All the grids produced in this work were done with Gambit, which is the Fluent, Inc. proprietary grid 
generator. There are other grid generators available with their own distinct advantages, but despite the bugs 
encountered in Gambit, we found it to be a useful tool. 

First, the domain is scaled up by 1000 to compensate for the 32-bit precision available in the grid generator. 
The Stirling engine, when entirely modeled, has a fairly wide range of length scales ranging from a 
thousandth of an inch in the narrow seal and appendix gap region to an inch or more in the heat 
exchanger regions. This scaling has the effect of making very small numbers larger so that smoother 
grids can be created. Once the mesh is created satisfactorly, the entire grid domain is scaled back down 
by a 1000 before utilizing it. 

Second, try to use structured or paver unstructured grids wherever possible to reduce the number of bad 
cells and to reduce the number of cells required. Avoid unstructured tetrahedrals because of the 
considerably reduced rate of solution convergence. 
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Third, be sure a sliding interface has essentially matching cell sizes on both sides of the interface to prevent 
interpolation errors and avoid instabilities. 

Fourth, be sure to place grid points in the boundary layer to properly model friction losses. 

Fifth, use sizing functions to smoothly expand the cell size from the boundaries at a rate of 10 percent. 

Sixth, use layering wherever possible, but avoid grid adaption and remeshing because this slows down the 
solution time considerably. Choose a layering size large enough to permit a reasonable 160 time steps 
per cycle, but small enough to allow for smooth transitioning to boundary layers. 

Seventh, utilize moving boundary layers on all walls that are moving to enhance conjugate heat transfer. 
And place a boundary layer on both sides of the interface (on the solid side and the fluid side). 

Eighth, include as much geometry as is practical to avoid having to redo the grid later when simplifying 
assumptions are found to matter. 

Ninth, clean up the geometry first in a CAD package to avoid creating virtual geometry that is not as 
robust. 

Tenth, test your grid by forcing the piston and displacer to be stationary and simulate the heat transfer in 
the engine. Once a converged solution is observed, continue running the solver and look for grid induced 
heat sources. It should be noted that Cleveland State University found producing a steady-state result 
first provided a better starting condition for the transient solution. But this last suggestion goes a bit 
further in providing a means for improving the grid by ” over- converging” the steady solution. 

For the successful case presented here, we used about 700,000 grid points and all grid boundary layers 
start at a distance of .001 in from the fluid/solid interface. 

In Fig. 11, we show the entire engine, piston region and spider region grid utilized in this simulation. 
Notice the piston has a moving boundary layer and grid layering is used to avoid small cells when the piston 
fully compresses. Also, notice in the spider region that all helium/wall interfaces have boundary layers on 
both sides to enhance conjugate heat transfer. The overall engine can be seen to not include any grid inside 
the heater or cooler since the heat exchanger wall temperture is specified as constant. Also notice there is 
some extra wasted grid in the spider region near the compression space. This is necessary to maintain grid 
quality in the other nearby boundary layer regions. Grid quality is more imporant than grid count for faster 
conjugate heat transfer calculations. 

In Fig. 12, the hot half of the engine is shown, including the gridding in the displacer and appendix gap. 
Notice that despite the complexity, tetrahedral grids were not required. The layering inside the displacer 
was done to avoid moving all the grid points inside the displacer for faster run-times. For example, without 
grid layering (the addition or subtraction of rows of cells) , all the grid cells inside the displacer would need to 
move. Instead, only the single layer of cells located at the top of the displacer (right side in picture) moves 
at a time in Fig. 12(b). Also, the expansion space uses layering to avoid skewed cells when the displacer is 
at top dead center (TDC). 

In Fig. 13, the gridding utilized in the seals and compression space are shown. The challenge in the seal 
region is to avoid overly skewed cells, overly mismatched cells across the interface, and reasonable cell size 
variation near the boundary layer. It is quite difficult to balance all these constraints. Layering was used 
in the compression region to avoid crushed cells and the axial direction spacing was chosen to maximize the 
time- step. 

In Fig. 14 the seal inside the displacer is shown. Notice the use of unstructured quadrilaterals throughout 
to avoid the inefficient tetrahedral gridding. Also notice just how small the seal gap is even compared to 
the thickness of the displacer wall. This region was very difficult to grid and required utilizing suggestion 10 
in Section C. Velocities are high in this region and this grid may not be adequate here, but since this seal 
region is not being studied this was left as seen. 
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(a) Entire Engine 


(b) Piston Region 



(c) Spider Space 


Figure 11. Grid Strategy 
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(a) Right Half 


(b) Layering Inside Dome 



(c) Layering in Expansion Space 


(d) Appendix Gap 



(e) Appendix Gap Close 


Figure 12. Hot End Grid Strategy 
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(a) Top and Bottom Seal, Rod, Compression Space 


(b) Top Seal and Compression 



(c) Top Seal Close-up 


(d) Bottom Seal and Compression 


Figure 13. Seal Grid Strategy 
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(a) Left Half 



(b) Inside Seal 



(c) Inside Seal Close-up 


(d) Inside Seal Close-up2 


Figure 14. 


Inner Seal Grid Strategy 
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D. Solver Optimization 


While one’s choices are limited with commercial codes, they do offer a variety of low-order solvers to choose 
from, including low Mach number variable density flows. 115 A comparison of two commercial codes with the 
latest high-order numerical techniques demonstrates some of the differences a solver decision can make. 50 
One of the limitations currently in Fluent is only first order accuracy is available in time when the mesh 
moves. CFD-ACE permits a higher accuracy in time with moving meshes, but does not have a coupled 
solver. The results shown in the next section utilize the segregrated implicit solver in Fluent, although it 
was later discovered the coupled implicit solver is potentially significantly better. Of course, the grid used for 
the segregrated solver may need further improvement if a coupled solver is used. It was found that the best 
method in Fluent and the best method in CFD-ACE performed comparably on simpler heat transfer test 
problems. 50 And with sufficient grid quality and density, the answers are comparable to the latest numerical 
methods. 

One of the challenges to utilizing actual engine geometry is the difficulty of creating a high quality mesh 
as discussed previously. However, certain switches in the solvers may be adjusted to attempt to compensate 
for poor grid quality in regions of conjugate heat transfer. 116 In what follows are specific recommendations 
for improving conjugate heat transfer and they contain Fluent specific switches. The other commercial 
packages should be modified with their own specific format. 

First, skewed meshes require secondary temperature gradients to be turned off to improve stability and 
to reduce undershoot /overshoot of local cell temperatures. In Fluent, this is accomplished by the 
command, ’’(rpsetvar ‘temperature/secondary-gradient? #f)” 

Second, in conjugate heat transfer problems with vastly different thermal conductivities between the fluid 
and solid zones, cells near the interface can experience round-off errors. To mitigate this effect, we 
want to ensure that the coarse levels of the AMG are visited. It is recommended that the default AMG 
cycle for energy be changed from Flexible to either F, V, or W-cycle. The F and V-cycles are more 
efficient in parallel simulations. The termination criteria should also be reduced to either 0.01 or 0.001 
in order to tighten up convergence. 

Third, although a reasonable temperature distribution can be achieved using single precision, in order to 
ensure that the total heat transfer is balanced the use of double precision is recommended. 

Fourth, if convergence still remains problematic with secondary gradients turned off, it is recommended 
that explicit under-relaxation of temperature be used to stabilize the solution. When this is activated, 
the energy equation is solved with a URF of 1, but the solution is explicitly relaxed with respect 
to the previous iteration. This allows for the equations to be efficiently solvered (URF equal 1) and 
instabilities due to poor quality mesh and non-linear material properties to be mitigated (relaxation 
factor). The GUI entry of URF for energy will specify the explicit relaxation factor to be used. To 
active this, the scheme command is: (rpsetvar Temperature/explicit-relax? T^t) 
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VII. Comparison with Analytical and Experimental Data 


The approach just described was applied to simulating a Stirling engine built by Infinia Corp. 117 (formerly 
Stirling Technology Company) being considered for use in space missions. The engine is currently being tested 
for long-duration missions and hence a good deal of experimental performance information is available for 
comparison. The geometry is export controlled for these engines and so limited information may be shared 
with the public. The entire simulation was performed first, with no knowledge of the ” correct” answer either 
from experiment or calibrated one-dimensional analysis. The operating conditions used in the simulation 
are shown below (Table 1):. 


Hot-End Temperature, K 

923 

Cold-End Temperature, K 

353 

Ambient Temperature, K 

293 

Frequency, Hz 

80 

Mean Pressure, Pa 

2.429E+6 

Power Piston Amplitude, mm 

6.0e-3 


Table 1. Operating Conditions 

The pressures and heat transfer occuring in various parts of the engine over two cycles are shown in 
Fig. 15(a). The pressure in the left (cooler end), middle, and right (heater end) of the regenerator is shown 
in Fig. 15(b). The average pressure in the expansion and compression spaces are shown to be about equal 
in Fig. 15(c). And, the temperature after 412.376 cycles of simulation are shown in Fig. 16 for an instant in 
time (135.36 degree Crank angle with reference to 0 degree displacer located at mean position) in the cold 
end of the machine and for the hot end are shown in Fig. 17. Notice the effects of radial heat transfer on 
the regenerator. 

An initial comparison with Sage suggested the axisymmetric simulation was over-predicting the indicated 
power as shown in Table 2. And the overall heat transfer in and out fell within a band predicted by Sage in 
which heat loss to the environment was included and then assumed excluded (adiabatic). 



Axisymmetric 

Simulation 

Modified Skupinski Sage ID 
With Ambient Heat Loss 

Modified Skupinski Sage ID 
No Ambient Heat Loss 

PV Power, W 

79.65 

68.77 

70.14 

Heat In, W 

247.315 

260.94 

193.8 

Heat Out, W 

168.313 

191.9 

123.7 

PV Efficiency 

.322 

.264 

0.362 

Pressure Ratio 

1.209 

1.187 

1.187 

Regenerator Delta Press., Pa 

10466.3 

15,888 

15790 

Heater Delta Press., Pa 

682 

194 

192 

Cooler Delta Press., Pa 

896 

316 

315 

Pressure Amplitude., Pa 

225269 

203100 

203200 

Mean Pressure, Pa 

2429130 

2378000 

2378000 


Table 2. Comparison With Sage ID Results 
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Absolute Pressure (Pa) 




(a) Pressure and Heat Transfer over Two Cycles 


(b) Pressure Along Regenerator over Two Cycles 



(c) Comparison of Expansion and Compression Pressure 


Figure 15. Pressure and Heat Transfer over Two Cycles 
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(a) Left Half 


(b) Compression and Cooler 



(c) Cooler and Seal 


(d) Seal Region 



Figure 16. 


Cold End Temperature 
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(a) Right Half 


(b) Below Heater 
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5.24e+02 
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4.56e+02 
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4.22e+02 
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3.53e+02 


WM 






Contours of Static Temperature (k) (Time=5.1 547e+00) 


(c) Around Heater 


(d) Regenerator 


Figure 17. 


Hot End and Regenerator Temperature 
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But then we gathered the actual performance of two identical engines identified as Technology Demon- 
strator Convertor (TDC) #13 and #14 from over many thousands of hours of operation since August 2003. 
As shown in Table 3, the axisymmetric solution which assumed no ambient heat loss to the environment, 
predicted not only the higher indicated power more accurately, but the efficiency was essentially within the 
variability of each built machine. 

A few words of caution about this encouraging result are in order. First, due to experimental limitations 
the piston amplitude is not known for certain. The original designs call for a piston amplitude of 6+/- .08 mm. 
And the mean operating pressure is 2.517 Mpa compared to the simulated mean pressure of 2.42913 Mpa. 
If the piston amplitude in the experiment is exactly 6.00mm, then this mean pressure difference implies the 
simulated power would need to be corrected up to 82.5874W. But since the engines are tuned for maximum 
performance and the energy measurements are considered more reliable than the piston amplitudes, it seems 
appropriate to conclude this simulation did in fact predict the performance and efficiency within the bounds 
of variability of the two engines tested. Since a slightly lower piston amplitude would counter-balance the 
corrected higher power just shown, it is important to improve the experimentally derived piston amplitude 
for future comparisons. 



T h 

T c 

freq 

net Qin 

PV Pout 

net Qout 

TDC #13 

646.0 

92.4 

81.4 

242.1 

78.2 

163.9 

TDC #14 

646.5 

94.4 

81.4 

250.4 

79.6 

170.8 

Simulation 

650.0 

80.0 

80.0 

247.3 

79.7 

168.3 

%Err-AVE 

100.6% 

100.3% 

98.2% 

100.4% 

100.9% 

100.6% 

%Err-TDC#13 

100.5% 

100.2% 

98.2% 

102.1% 

101.9% 

102.7% 

%Err-TDC#14 

100.6% 

100.3% 

98.2% 

98.8% 

100.0% 

98.5% 


Table 3. Comparison With Experimental Results 


VIII. Conclusion 

For the first time in the U.S., a fully converged axisymmetric simulation of an actual Stirling engine was 
successfully performed. And moreover, by utilizing the latest advances in computer processors, communi- 
cation hardware, parallel decomposition strategy, dynamic meshing, solution algorithms, and careful grid 
design, it is possible to compute one complete Stirling cycle in less than an hour. 

While not conclusive until additional testing is completed, the full axisymmetric simulation predicted 
the engine performance within the bounds of both engine test’s data. And predicted the overall efficiency 
and indicated power considerably better than our initial efforts with using a leading one-dimensional Stirling 
analysis prediction on the first attempt. The one-dimensional analysis can certainly be improved by mod- 
ifying various friction and heat transfer factors, either from experimental data or from this axisymmetric 
solution. But the real value of the axisymmetric solution is no adjustments were required because the engine 
is simulated from first principles. Of course the one-dimensional analysis is valuable for other reasons such 
as fast optimizations and performance mapping. 

These results and observations lead to the conclusion that one, two, and three-dimensional modeling 
should all be employed and all three paradigms provide important capabilities that when combined provide 
potent combination of initial design, empirical coefficient adjustment, optimization, and final prototype 
demonstration before the first metal is cut. 
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